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Abstract 

We give a systematic method for discretizing Hamiltonian partial differential 
equations (PDEs) with constant symplectic structure, while preserving their 
energy exactly. The same method, applied to PDEs with constant dissipa- 
tive structure, also preserves the correct monotonic decrease of energy. The 
method is illustrated by many examples. In the Hamiltonian case these in- 
clude: the sine-Gordon, Korteweg-de Vries, nonlinear Schrodinger, (linear) 
time-dependent Schrodinger, and Maxwell equations. In the dissipative case 
the examples are: the Allen-Cahn, Cahn-Hilliard, Ginzburg-Landau, and heat 
equations. 

Keywords: Average vector field method, Hamiltonian PDEs, dissipative 

PDEs, time integration. 
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1. Introduction 

"The opening line of Anna Karenina, 'All happy families resemble one 
another, but each unhappy family is unhappy in its own way', is a useful 
metaphor for the relationship between computational ordinary differential equa- 
tions (ODEs) and computational partial differential equations (PDEs). ODEs 
are a happy family - perhaps they do not resemble each other, but, at the very 
least, we can treat them by a relatively small compendium of computational 
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techniques. . . PDEs are a huge and motley collection of problems, each unhappy 
in its own way" (Quote from A. Iserles' book |15j). 

Whereas there is much truth in the above quote, in this paper we set out to 
convince the reader that, as far as conservation or dissipation of energy is con- 
cerned, many PDEs form part of one big happy family (cf. also [17 ) that, after 
a very straightforward and uniform semi-discretization, may actually be solved 
by a single unique geometric integration method - the so-called average vector 
field method - while preserving the correct conservation, respectively, dissipa- 
tion of energy. The concept of 'energy' has far-reaching importance throughout 
the physical sciences [TO]. Therefore a single procedure, as presented here, that 
correctly conserves, resp. dissipates, energy for linear as well as nonlinear, low- 
order as well as high-order, PDEs would seem to be worth while. 

Energy-preserving schemes have a long history, going back to Courant, Fricdrichs, 
and Lewy's cunning derivation [6] of a discrete energy conservation law for 
the 5-point finite difference approximation of the wave equation which they 
used to prove the scheme's convergence. The conservation law structure of 
many PDEs is considered fundamental to their derivation, their behaviour, and 
their discretization. Li and Vu-Quoc [H] give a historical survey of energy- 
preserving methods for PDEs and their applications, especially to nonlinear 
stability. What is relevant to us here is that many of these methods (e.g. 
[H [H HH Uni HOI [HI EH! ESI) have an ad-hoc character and are not com- 
pletely systematic either in their derivation or in their applicability; in contrast, 



the method discussed here (Eq. (16) below) is completely systematic, applies 
to a huge class of conservation and dissipative PDEs, and depends functionally 
only on the PDE itself, not its energy. In some cases it reduces to previously 
studied methods, for example, it reproduces one of Li and Vu-Quoc's schemes 
[T5| for the nonlinear wave equation. Even in these cases, however, it sheds 
considerable light on the actual structure of the scheme and the origin of its 
conservative properties. See also the discussion comparing different construc- 
tions of energy-preserving integrators in [7]. 

We consider evolutionary PDEs with independent variables (x, t) £ M. d x K, 
functions u belonging to a Hilbert space B with valuefQu(x, t) £ K m , and PDEs 
of the form 

(1) 

where T> is a constant linear differential operator, the dot denotes J^, and 

H[u] = [ H{x-u {n) ) dx (2) 
Jn 

where f2 is a subset of R d x R, and dx = dx\dxi . . . dxa- is the variational 



1 Although it is generally real- valued, the function u may also be complex- valued, for ex- 
ample, the nonlinear Schrodinger equation. 
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derivative of % in the sense that 

d 



de 



H\u- 



I e=0 



sn 

Su ' 



(3) 



for all u,v £ B (cf. [21]), where (,) is the inner product in B. For example, if 
d = m = l,B = L 2 (ft), and 



H[u]= / H(x;u,u x ,u xx ,...)dx, 



then 



Sn_cUJ_ d fdH\ 2 / dH 
Su du x \du x J x \du xx 



(4) 
(5) 



when the boundary terms are zero. 

Similarly, for general d and to, we obtain 



sn 

Sui 



dH 

dui 



E 

it=i 



d ( dH 



dx k \dui t k 



1 = 1, 



(6) 



We consider Hamiltonian systems of the form Q, where I? is a constant 
skew symmetric operator (cf. [28]) and T-L the energy (Hamiltonian). In this 
case, we prefer to designate the differential operator in ([lj with S instead of T>. 
The PDE preserves the energy because S is skew-adjoint with respect to the L 2 
inner product, i.e. 

(7) 



iSudx = 0, Vtt G B. 



The system \1\ has 1 : B 



I as an integral if I = L P<S^ dx = 0. 
Integrals C with T>¥^- — arc called Casimirs. 

° ou 

Besides PDEs of type (IT]) where T> is skew-adjoint, we also consider PDEs 
of type (nj) where I? is a constant negative (semi)dcfinite operator with respect 
to the L inner product, i.e. 



uDu dx < 0, Vu € B. 



(8) 



In this case, we prefer to designate the differential operator V with AT and the 
function T-L is a Lyapunov function, since then the system ([TJ), i.e. 



sn 

Su ' 

Jfl OU Oil 



(9) 



has n as a Lyapunov function, i.e. n = j n j^Nj^ dx < 0. We will refer 
to systems (JTJ) with a skew-adjoint S and an energy n as conservative and to 
systems (JTJ) with a negative (semi) definite operator M and a Lyapunov function 
n as dissipative. Note that the operator N need not be self-adjoint. (In Example 
10 the Ginzburg-Landau equation, Af = d x + ed x is not self-adjoint.) 
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Conservative PDEs ([I]) can be semi-discretized in "skew-gradient" form 

u = SVH{u), S T = -S, ueR k (10) 

when T> = S is skew-adjoint. Here, and in the following, we will always denote 
the discretizations with bars. H : R k — s- R is chosen in such a way that HAx is 
an approximation to H. 

Lemma 1. Let 

H[u] = I H(x;u {n) )dx, (11) 
Jn 

and let HAx be any consistent (finite difference) approximation to H (where 
Ax := AxiA:c2 . . . Axd) with N degrees of freedom. Then in the finite- dimensional 
Hilbert space M. N with the Euclidean inner product, the variational derivative 

— (HAx) is given by VH. 
ou 



Proof. We denote the consistent (finite difference) discretization of H[u] = 

J n H(x; «(")) dx by H (u) — H(u)Ax where u £ R N denotes the discrete values 
of u, in the multidimensional case after choosing an ordering. The variational 

, • • m ■ , 

derivative — — is then given by 

ou 




It holds for the directional derivative that 

d 



de 



H{u + ev)\ e=Q = (v-VH). (13) 



Since Eqs ( 12 ) and ( 13 ) must hold for all vectors v, we have 



AT-/ 

— Ax = V7Z', and hence — = VH. (14) 
ou ou 

It is worth noting that the above lemma also applies directly when the ap- 
proximation to H is obtained by a spectral discretization, since such an ap- 
proximation can be viewed as a finite difference approximation where the finite 
difference stencil has the same number of entries as the number of grid points 
on which it is defined. 

The operator V is the standard gradient, which replaces the variational 
derivative because we are now working in a finite (although large) number of 
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dimensions (cf. e.g. (|6|)). When dealing with (semi-) discrete systems we use 
the notation Uj t7l where the index j corresponds to increments in space and 
n to increments in time. That is, the point u^ n is the discrete equivalent of 
u(a + jAx, t + nAt) where x g [a, b] and where to is the initial time. In most 
of the equations we present, one of the indices is held constant, in which case, 
for simplicity, we drop it from the notation. For example, we use Uj to refer to 
the values of u at different points in space and at a fixed time level. 

Theorem 1. Let S (resp. M) be any consistent constant skew (resp. negative- 
definite) matrix approximation to S (resp. M). Let HAx be any consistent 
(finite difference) approximation to %. Finally, let 

f(u):=SVH(u) {resp. f(u) :=AfVH{u)), (15) 

and let u n be the solution of the average vector field (AVF) method 



l n+l 



At 



f + (16) 

Jo 



applied to equation (15). Then the semidiscrete energy H is preserved exactly 



(resp. dissipated monotonically): 

H(u n+ i) = H(u n ) (resp. H(u n+1 ) < H(u n )). 
V. is preserved by the flow of ii = <SVH(m) since 

H = (VH) T SVH = 0. (17) 

Discretisations of this type can be given for pseudospectral, finite-element, 
Galerkin and finite-difference methods (cf. [24l [25]); for simplicity's sake, we 
will concentrate on finite-difference methods, though we include two examples of 
pseudospectral methods for good measure. We consider only uniform grids; see 
|16j for finite difference discretizations of this type on nonuniform grids, which 
require inner products (u n ,v n ) = u n ■ (Mv n ) with M =/= I. 

The AVF method was recently [53] shown to preserve the energy H exactly 
for any vector field / of the form f(u) = SVH(u), where % is an arbitrary 
function, and S is any constant skew matrbj^] The AVF method is related 
to discrete gradient methods (cf. [53]); amongst discrete gradient methods it 
is distinguished by its features of linear covariance, automatic preservation of 
linear symmetries, reversibility with respect to linear reversing symmetries, and 
often by its simplicity. It is one member of the family of Galerkin methods 
introduced by Betsch and Steinmann [2] for the case of canonical classical me- 
chanical systems, although the dependence of the method on / alone (and not 
%) was not realized there. 



2 The relationship of (1 1 6 h to Rungc-Kutta methods was explored in [5]. 
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We remark that although ( 16 ) does not specify the method in a completely 
closed form (because of the integral), the same is true for implicit methods, 
whose description needs to be supplemented by an iterative solver in any im- 
plementation. Despite this, broadening the class of methods from explicit to 
implicit is essential for many differential equations and allows new properties 
such as A-stability and symplecticity. Similarly, a further broadening to include 
integrals of the vector field allows the new properties of energy conservation and 
dissipation. 

In many cases — for all linear, polynomial, or scalar terms in /, which includes 
all the examples in this paper — the integral can be evaluated exactly in closed 
form. In other cases, it can be approximated by quadrature to any desired 
degree of accuracy. 

If V is a constant negative-definite operator, then the dissipative PDE ([I]) 
can be discretized in the form 

u =JfVH{u), (18) 

where Af is a negative (semi)definite matrix and % is a discretization as above. 
That is, W is a Lyapunov-function for the semi-discretized system, since 

jl= (VH) T 7VW< 0. (19) 



The AVF method ( 16 ) again preserves this structure, i.e. we have 

n(u n+1 ) <n(u n ), (20) 

an d % is a Lyapunov function for the discrete system. Taking the scalar product 
of (16 1 with J Q VH((1 — u n + £ u n+i) d£ on both sides of the equation yields 



1 

At 



i.e. 



/ {U n+X - Un) • VH((1 - 0«n + £«n+l) < ( 21 ) 

Jo 

/ — U((l-0u n + tu n+1 )dl;<0, (22) 
Jo "? 



j_ r 1 d 

At 



and therefore 



^(% +1 )-%))<0. (23) 

Our purpose is to show that the procedure described above, namely 

1. Discretize the energy functional H using any (consistent) approximation 
HAx 

2. Discretize I? by a constant skew-symmetric (resp. negative (semi) definite) 
matrix 

3. Apply the AVF method 



ENERGY CONSERVING PDES 



7 



can be generally applied and leads, in a systematic way, to energy-preserving 
methods for conservative PDEs and energy-dissipating methods for dissipative 
PDEs. We shall demonstrate the procedure by going through several well-known 
nonlinear and linear PDEs step by step. In particular we give examples of how 
to discretize nonlinear conservative PDEs (in subsection 2.1), linear conservative 
PDEs (in subsection 2.2), nonlinear dissipative PDEs (in subsection 3.1), and 
linear dissipative PDEs (in subsection 3.2). 

Energy dissipation has been much less studied than energy conservation. 
Stuart and Humphries |33j consider conditions for certain linear methods (such 
as backward Euler) to be dissipative for gradient systems x — — VF(x), F(x) > 

0, lim|| ^ || j.Q F(x) = oo. This class is however far more restrictive than the 

dissipative systems the restriction from Af negative (semi) definite and not 
necessarily self-adjoint to Af — —id alone is substantial [23] ■ In particular, 
backward Euler is not dissipative for systems of the form Q (see Example 10 
below) . 



The method ( 16 ) is formally second order in time. The relationship between 
accuracy in space and time is a complicated one that we cannot explore here. 
Depending on the PDE and the scientific goals the accuracy in time may need to 
be much less, the same, or much greater than that in space. The temporal order 
can be increased if necessary by composition |22| or by including derivatives of 
the right hand side [2~§j . 

Steps 1 and 2 yield semidiscretizations that are Hamiltonian (or Poisson) in 
the conservative case and dissipative in the dissipative case. Thus they could be 
followed by a symplectic time integrator, so this procedure unifies symplectic and 
energy-conserving integration of conservative PDEs and allows for systematic 
comparisons between the two. Other approaches to energy-conserving integra- 
tion, e.g. [IH1 HOI H5 , conflate time and space discretization and obscure the 
fundamental role played by the semidiscrete energy and its Euclidean gradient. 

The choice of a symplectic versus an energy-preserving time integrator is an 
interesting and delicate question that depends on the PDE and the scientific 
goals. Some factors favour energy conservation: 

1. When the (semidiscretized) energy level sets are compact, conserving en- 
ergy may (via energy stability) improve performance at large time steps; 
Simo and Gonzalez [32 show that the unresolved high frequencies are con- 
trolled by exact energy conservation, whereas they can lead to instability 
in the symplectic midpoint rule. (See also their comparisons in the context 
of relative equilibria [12] and in nonlinear elastodynamics |31j.) 

2. If a nonconservative scheme is used for a system of conservation laws and 
large (e.g.) mass errors result, this is not taken to be a sign that the spa- 
tial mesh size should be reduced, or that mass should be globally rescaled 
back to its original value; rather, it is taken to be a sign that a conser- 
vative scheme should have been used. Unitarity in Schrodinger equations 
and orthogonality in spin chains followed a similar course: their preser- 
vation was first discovered to be important; then enforced by projection; 
eventually, intrinsically-conserving schemes were discovered that are now 
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widely used. 

3. In fully chaotic dynamics, symplecticity may be irrelevant. It is not used in 
statistical mechanics (only energy and volume preservation are required). 
If the flow on an energy level set is Anosov, then it is structurally stable, 
even to non-Hamiltonian perturbations. 

4. It is easier to adapt the time step in energy-preserving than in symplcctic 
integration. 

5. It leads to different constraints on the eigenvalues and bifurcations of or- 
bits than symplecticity; for example, energy preservation leads to periodic 
orbits occurring in 1-parameter families (parameterized by the energy), 
whereas symplecticity does not. 

Some factors favour symplecticity: 

1. Unlike energy, symplecticity is the defining property of Hamiltonian me- 
chanics. It ensures preservation of many phase space features such as 
genericity of quasiperiodic orbits. 

2. It can be essential to capture qualitatively correct long-time limit sets and 
long-time statistics. 

3. In 2n dimensions, symplecticity provides n(2n — 1) constraints, energy 
only one. 

4. It gives the user a very convenient check on the simulation, namely, the 
energy error. However, the security this provides may be illusory as neither 
conservation of energy nor small energy errors in a nonconservative scheme 
ensure that the simulation is reliable. 



2. Conservative PDEs 



2.1. Nonlinear conservative PDEs 

Example 1. Sine-Gordon equation: 
Continuous: 

d 2 tp d 2 ip 

The Sine-Gordon equation is of type ([I]) with 



dip 
Ox 



+ a (1 — cos</?) 



dx, 



(24) 



(25) 



where u :- 



<P 



and 



S 



1 

-1 



(26) 



(Note that it follows that n = ^.) 

Boundary conditions: periodic, u(—20,t) — u(20,t). 



ENERGY CONSERVING PDES 



9 



Semi-discrete: finite difference^ 



2(Ax) 



S 



id 
-id 



The resulting system of ordinary differential equations is 



7T 



= SVH fd = 
where L is the circulant matrix 



^2 Lip -asinip 



(27) 
(28) 

(29) 



L 



-2 1 



1 



We have used the bold variables cp and ir for the finite dimensional vectors 
[ipi, (f2, ■ ■ ■ , </5jv] T , et cetera, which replace the functions tt and ip in the (semi-) 
discrete case. Where necessary, we will write ip n , et cetera to denote the vector 
cp at time to + nAt. 

The integral in the AVF method can be calculated exactly to giv^] 



I 

At 



7T,, 



<Pr. 
7Tk 



(30) 



("Tr, 



w n )/2 



i (V„+i+V„)/2-a(cos ^ 



cos Vn)/Or, 



Semi-discrete: spectral discretization 

Instead of using finite differences for the discretization of the spatial deriva- 



tive in (25), one may use a spectral discretization. This can be thought of as 



replacing <p with its Fourier series, truncated after N terms, where N is the 
number of spatial intervals, and differentiatingthe Fourier series. This can be 
calculated, using the discrete Fourier transform 5 ! (DFT) , as T^ 1 DnFn<P where 
J-n is the matrix of DFT coefficients with entries given by [J-N]n,k = uT 



nk 
N ) 



3 Summations of the form mean J^jLo 1 umess stated otherwise. 

4 For numerical computations, care must be taken to avoid problems when the difference 
<p n +l — <p n in the denominator of (|30 b becomes small. We used the sum-to-product identity 
cos a — cos b = — 2 sin((a + b)/2) sin((a — b) /2) to give a more numerically amenable expression. 
In Eq. | |30| , elementwise division of vectors is used. 

5 In practice, one uses the fast Fourier transform algorithm to calculate the DFTs in 
O (TV log TV) operations. 
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lun = e l27T / N Additionally, [^jy 1 ]™,*: = jj ljJ N nk an d -Djv is a diagonal matrix 
whose (non-zero) entries are the scaled wave-numbers^] 



diag(£> 



N) 



2ni 



0,1,2,. 



N - 1 N- 1 

2 ' 2~ : 



-2,-1 



(for TV odd), where / = 6 — a is the extent of the spatial domain; that is l/N = 
Ax. (For more details on properties of the DFT and its application to spectral 
methods see [3] and [53].) 



\A + \ [ :f n 1d nFn<p\] + a(l - cos^-) 



5 



id 

-id 



The resulting system of ODEs is then given by 



<P 

7T 



— SVn 8 j, — 



7T 



Again, the integral in the AVF method can be calculated exactly to give 

f n+l - <Pr, 



At 
At 



(TTn+l + 1T„)/2, 



(31) 
(32) 

(33) 
(34) 



- a(cosy>, i+1 - cos (p n )/(<f n+1 - V„)- ( 35 ) 

Initial conditions and numerical data for both discretizations: 
Spatial domain, number N of spatial intervals, and time-step size At used 
wer^H 



xe [-20,20], N = 200, 
Initial conditions: 

<p{x,0) = 0, 

7r(x,0) = 



cosh(2ir) 



At = 0.01, parameter: a = 1. 



Right-moving kink 

and left-moving (36) 

anti-kink solution. 



6 Care must be taken with the ordering of the wave numbers since different computer 
packages use different effective orderings of the DFT/IDFT matrices in their algorithms. 
Additionally, one must ensure that all modes of the Fourier spectrum are treated symmetrically 



— for N even, this requires replacing the k - 
7 Here and below, if x S [a, b], then Aa; = 



entry with zero to give [0, . . . , 
, and Xj = a + j Ax, j = 0, 1, 



N . 1Q -N 



, N. 
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— Midpoint 
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time 
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SG eqn, Nspace= 201, At = 0.01 




Figure 1: Sine-Gordon equation with finite differences semi-discretization: Energy error (left) 
and global error (right) vs time, for AVF and implicit midpoint integrators. 




Figure 2: Sine-Gordon equation with spectral semi-discretization: Energy error (left) and 
global error (right) vs time, for AVF and implicit midpoint integrators. 



Numerical comparisons of the AVF method with the well known (symplec- 
tic) implicit midpoint integratoij^] are given in figure [l] for the finite differences 
discretization, and in figure [2] for the spectral discretization. 

Example 2. Korteweg-de Vries equation: 



Continuous: 



du du d 3 u 

dt dx dx 3 ' 



H = 



(u x ) - u 3 



dx. 



(37) 
(38) 



'Recall that the implicit midpoint integrator is given by 
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S 



dx 



(39) 



Boundary conditions: periodic, u(— 20, t) = u(20,t). 
Semi-discrete: 



i 



2{Ax)'- 



I \ 2 3 



(40) 



S = 



2Ax 



-1 

1 -1 



(41) 



1 -1 
1 



Initial conditions and numerical data: 



x e [-20, 20], TV = 400, At = 0.001. 

Initial condition: u(x,0) = 6 sech 2 (a:) (for two solitons). Numerical compar- 
isons of the AVF and the midpoint rule are give in Figure [3] The global errors 
are comparable. The AVF is particularly easy to implement here because the 



integral ( 16 ) is evaluated exactly by Simpson's rule. 



Example 3. Nonlinear Schrddinger equation: 
Continuous: 



d 



i 
-i 



dt V W 

where u* denotes the complex conjugate of u 



n = 



m 1 1 

8u* 



1 




du 


2 + ?m; 






dx 



dx, 



S = 



i 
-i 



Boundary conditions: periodic, u(— 20, £) = it(20, t). 
Semi-discrete: 



i 



(Ax] 

5 = i ' -id 
Initial conditions and numerical data: 



2 - u 

id 



2 . 7, ,4 

I 



(42) 

(43) 
(44) 



(45) 
(46) 
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Figure 3: Korteweg-de Vries equation: Energy error (left) and global error (right) vs time, 
for AVF and implicit midpoint integrators. 
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Figure 4: Nonlinear Schrodinger equation: Energy and global error vs time, for AVF and 
implicit midpoint integrators. 
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10 -s NLS equation, Nspace= 201, At = 0.05 




Figure 5: Nonlinear Schrodinger equation: Total probability error vs time, for AVF and 
implicit midpoint integrators. 



x G [-20,20], 
Initial conditions: 



N = 200, 



$lu(x,0) 
^su(x, 0) 



At = 0.05, 



parameter: 7 = 1. 



exp (-Or - l) 2 /2) 
exp (-x 2 /2) . 



(47) 



Numerical comparisons of the AVF and the midpoint rule are give in Figur es |4} 
[5j The global errors are comparable. As for the KdV equation, the integral ( |16| ) 
is evaluated exactly by Simpson's rule. The AVF method does not conserve the 
total probalility because it is not unitary. 

Example 4. Nonlinear Wave Equation: 

Continuous: 

The 2D wave equation 

^ = Atp- 91 ^, <p = (p{x,y,t), (x,2/)e[-l,l]x[-l,l],t>0, (48) 
is a Hamiltonian PDE with Hamiltonian function 



H 



dx dy, 



(49) 



where 7r = d/dtip and the operator S is the canonical 2x2 symplectic matrix. 
Boundary conditions: periodic. 
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Semi-discrete: spectral elements 

We discretize the Hamiltonian in space with a tensor product Lagrange 
quadrature formula based on p + 1 Gauss-Lobatto-Legendre (GLL) quadrature 
nodes in each space direction. We obtain 



, p p 



3 1=0 32=0 



where d 



3l, k 



dl k (x) 



dx 



f P \ 2 / P \ 2 1 \ 

X) d ii,fcVfc,i2 ) +( d h,m<Pji,m\ +2^31,32 

\k=0 / \m=0 / J 

(50) 

and lk{x) is the fc-th Lagrange basis function based 



on the GLL quadrature nodes xq, . . . ,x p , and with wq, . . . , w p the corresponding 
quadrature weights. The numerical approximation is 



(51) 



k=0 m=0 



and has the property ip p (xj 1 ,yj 2 ,t) = <+Pj 1 ,j 2 {t), so that the data can be stored 
in the (p + 1) X (p + 1) matrix with entries <fij lt j 2 - 
Initial conditions and numerical data: 



(x,y) g [-1.1] 2 , V(<p) 



9 



Initial condition: </?(a;,y,0) = sech(10x)sech(10y), ir(x, y,0) = 0. 

In figure [6] we show some snapshots of the solution. The energy error is 
shown in figure [7J 

2.2. Linear conservative PDEs 

Example 5. (Linear) Time- dependent Schrddinger Equation: 
Continuous: 



du d 2 u 



(52) 



This equation is bi-Hamiltonian, i.e. it has 2 independent symplectic struc- 
tures. The first Hamiltonian formulation has 



H x = 



p7T 




du 


1 — 7T 




dx 



- V(x) \u\ 



dx 



and 



Si = 



i 
-i 



The second Hamiltonian formulation has 



U-2 



\u\ 2 dx 



(53) 
(54) 
(55) 
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T = 



T = 3.1250 




T = 6.875 



T = 10 




Figure 6: Snapshots of the solution of the 2D wave equation at different times. AVF method 
with step-size At = 0.6250. Space discretization with 6 Gauss Lobatto nodes in each space 
direction. Numerical solution interpolated on a equidistant grid of 21 nodes in each space 
direction. 



0.5 


-0.5 

g " 1 
>» -1 "5 



x iq-11 2D Wave equation avf and ode15s 



0) 



-2 
-2.5 

-3 



10 



time 



Figure 7: The 2D wave equation {481 . MATLAB routine odel5s with absolute and relative 
tolerance 10~ 14 (dashed line), and AVF method with step size At = 10/(2 5 ) (solid line). 
Energy error versus time. Time interval [0,10]. Space discretization with 6 Gauss Lobatto 
nodes in each space direction. 
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and 



dl- V(x) 

-dl + V(x) 

Boundary conditions: periodic, u(— ir,t) — u(n,t). 
Scmi-discrctc: 



i 



Si = i 

The second semi-discretization is 

n 2 = J2 



\u j+1 -Uj\ 2 - V(xj)\uj\ 



id 
-id 



S 2 = i 



A 

-A 



where 



/ -2-V 1 
1 -2-V 1 



V 



1 \ 



1 



(56) 



(57) 
(58) 

(59) 
(60) 

(61) 



L ... 1 -2-V J 

Both discretizations result in the same semi-discrete system and the AVF method 
(which in the linear case coincides with the midpoint rule) therefore preserves 
both Hi and "H 2 , as well as the two symplectic structures. 

Initial conditions and numerical data: 

xe[-7r,7r], iV = 50, At = 0.1, V(x) = 1 - cos(a;). 
Initial conditions: 



$tu{x,0) = e~ ( i )2 , Qu(x,0)=0. 



(62) 



Example 6. Maxwell's Equations (Id): 

We first look at the one-dimensional Maxwell equation 
Continuous: 



d_ 
dt 



E 
D 







'dx 



'dx 




5H 
SE 
5H 
SB 



H(E 



E 2 + B 2 ) dx, 



(63) 
(64) 
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8 10 



Figure 8: Linear Schrodinger equation: Error in energy Hi Ax vs time, AVF method 



and 



s 



w- 

dx 

#- o 

ox 



(65) 



S is skew-adjoint on {(E, B) £ C 1 : E(0) = E(l) = 0} (and therefore on the 
Sobolev space Hq), 

Boundary conditions: 



E(0,t) = E(l,t) = 0, 
H(0,t) = |f(l,t) = 0. 



(66) 



Semi-discrete: 

We now obtain % by discretizing % in a simple way by applying the trape- 
zoidal rule to the integral (64) at the points Xj — j^j and dividing by Ax, that 
is 

( 1 \ 1 N ~ l ( 1 \ 1 
H(E 1 ,---E N _ 1 ,B ,--- ,B N )= J2 +C I B o + E c o^ +C 4^' 

(67) 

where we have already used that E(xo,t) = E(xN,t) — 0. The differential 
operator S is discretized with central differences yielding 



S 



2V-1,JV-1 q 



(68) 
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o 

UJ 



CD 
1= 
UJ 




0.4 0.6 
time 



Figure 9: One-dimensional Maxwell equation: energy error vs time, AVF integrator. 



where the (N — 1) x (N + 1) matrix G is given by 



-2 1 
-1 



G = 



2Aa 



-1 1 
-10 2 



(69) 



Initial conditions and numerical data: 

Note that the Neumann boundary conditions are satisfied at least to order 
1 in space, despite the fact that we only intended to somehow approximate the 
energy %. The numerical experiments confirm that the discrete energy HAx is 
preserved to machine precision. Figure [9] shows the error of the AVF method 
for the Maxwell equation with N = 100, At = 0.001, c — 1, and initial value 



E(x,0) = e- 100 (^) 2 , B(x,0) = e - 100 (*-5) 2 



Example 7. Maxwell's Equations (3D): 



Continuous: 

We consider Maxwell's equations in CGS units for the electromagnetic field 
in a vacuum 

— I" B 1 ° ~ CVX 1 I" B 1 f70^ 

dt E cVx E ^ ' 
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with the operator 



Vx 



dz 



'hi 



d_ 

dz 



-sr 

OX 



_d_ d_ 

dy dx 



(71) 



This equation has two Hamiltonian formulations of type ([!]). The first Hamil- 
tonian formulation has the helicity Hamiltonian 



Ux = J [c\ bT (VxB) + c^E T (V x E)j dxdydz 
(cf. [1]) and the operator 



Si 



-h 

h o 



(72) 



(73) 



where ^3 designates the 3x3 unit matrix. The second Hamiltonian formulation 
has the Hamiltonian 



n 2 = J [c^B T B + c^E T E^j dxdydz 
(cf. [14]) and the operator 

S 2 





Vx 



-Vx 




(74) 



(75) 



Boundary condition: periodic on the unit cube Q. 
Semi-discrete: 

On a regular grid with lexicographical ordering in every component of E 
(resp. B) and concatenating the discretized components gives one discrete 
vector Eh (resp. Bh), the operator Vx is represented by a matrix A. The 
discretization in the first case is given by 



Ux = c l -ElAE h + c l -BlAB h 
and the obvious discretization of S\. For the quadratic Hamiltonian, 



1 



1 



H 2 = c-EiE h +c-B 1 h B h 



and 



S 2 



—A 
A 



Both discretizations result in the same semi-discrete system 



" B h ' 




' -A' 




' B h ' 


. Eh _ 




A 







(76) 

(77) 
(78) 

(79) 
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'cs 




Figure 10: Three-dimensional Maxwell equation: plots of energies HiAx (dash-dot) and 
H2AX (dash) vs time. 



and the average vector field method preserves both Hi and "H 2 - 
Initial conditions and numerical data: 

Preservation of Hi and Hi is numerically confirmed by an experiment with 
random initial data on a regular grid with 30 points in every direction and 
constant c = 1. The result of the AVF method with At = 0.01 can be seen in 
figure [lO] 



3. Dissipative PDEs 

3.1. Nonlinear dissipative PDEs 
Example 8. Allen-Cahn equation: 



Continuous: 



du 

at 



= du xx + u — u 3 , d > 0, 



H 



1 j/ \ 2 1 2 1 1 4 

— a 11, u H — u 

2 y ' 2 4 



-1. 



Boundary conditions: Neumann, it x (0, t) = 0,^(1,^ = 0. 
Semi-discrete: 



(80) 

(81) 
(82) 



JV r 
j=0 



2(Ax) 



2 v"j+i 



u 3 -) 



1 2,1 

— u„- H — lt w 

2 4 



(83) 
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Af = -id. 
Initial conditions and numerical data: 



(84) 



x £ [0, 1], N = 100, At = 0.001, parameter: d = 0.001. 
Initial condition: u(x, 0) = cos(-7rx). 
Example 9. Cahn-Hilliard equation: 
Continuous: 



du d 2 
~dt = dx^ 



(pu + ru 3 + qu xx ) , 



n = 



-pu + -ru - -q{u x ) 

N = dl. 



dx, 



Boundary condition: periodic, u(Q,t) = u(l,t). 
Semi-discrete: 



1 2 , 1 4 1 9 / ^2 
_p„. +_ rUj . ____(„. +1 _„.) 



AT 1 



(Ax) 2 



-2 1 
1 -2 1 



1 



Initial conditions and numerical data: 



1 -2 1 

1 -2 



(85) 

(86) 
(87) 

(88) 

(89) 



i€[0,l], TV = 50, At = 1/1200, parameters:p= -1, q= -0.001, r= 1. 
Initial condition: 

u(sg,0) = 0.1 sin(27rx) + 0.01 cos(4ttx) + 0.06 sin(47T2:) + 0.02 cos(IOttx). 
Example 10. Ginzburg-Landau equation: 
Continuous: 

A Ginzburg-Landau equation arising in a model of traffic flow is given by 
du 



dt 



= (d x + edl) [6« + d 2 x u - u 3 } , e > 0, 



(90) 



and is a slight modification of the model considered in and [27j . The equa- 
tion can be written as 

m= M lu-> (91) 
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0.5 1 1.5 

time 



CH equation, Nspace= 51, At = 0.00083333 




-0.14 1 1 1 1 

0.5 1 1.5 

time 



Figure 12: Cahn-Hilliard equation: Global error (left) and energy (right) vs time, AVF 
integrator. 
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with 
and 



N = d x + edl 



n 



3u z 



1 fduV 1 



2 \dx 



4 



dx. 



Note that M is not self-adjoint. 

Boundary condition: u(±5,t) = and u xx (±5,t) = 0. 
Semi-discrete: 



7V-1 



3w- 



and 
where 



1 / ltj' + l — 

Ax 



7V = A + eB, 

1 

-1 1 



1 4 

4 3 



Un = 0. 



.4 



2Aa 



-1 1 

-1 



is a discretization of d x , and 



(Ax) 



-2 1 
1 -2 1 



1 -2 1 
1 -2 



(92) 
(93) 



(94) 
(95) 

(96) 



(97) 



is a discretization of d\. The average vector field method preserves the decay 
of function % in contrast to some standard integrators. 
Initial conditions and numerical data: 



x G [-5, 5], N = 50, At = 0.001, parameter: e = 0.001. 
Initial condition: u(x, 0) = fT 100 ^) 2 . 

In figure [i~3| w e compare the AVF method with the backward Euler method. 
For the latter, % appears to be monotonic increasing instead of monotonic de- 
creasing. This quickly leads to a qualitatively incorrect solution. Since backward 
Euler is one of the most famous dissipative methods for M = —id, this clearly 
shows that the AVF method is not a trivial extension. 
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time time 



Figure 13: Ginzburg-Landau equation: Plots of the energy function HAx computed with 
the AVF method (left) and backward Eulcr (right). Note that Backward Euler exhibits an 
increase in the energy instead of the correct decrease. 



3.2. Linear dissipative PDEs 
Example 11. Heat equation: 

Continuous: 



The heat equation 



du 

di 



is a dissipative PDE and can be written in the form ([!]), i.e. 



du _ 8%i 
dt 1 Su 



du _ ^ 8I-L2 
dt 2 5u 



(98) 



(99) 



with the Lyapunov functions Hi(u) = J Q \u\dx and Hiiu) = \u 2 dx and 
the operators N\ = — 1 and A/2 = d%, respectively. 

Boundary conditions: u(0, t) = u(l, t) = 0. 

Semi-discrete: 



and 



as well as 



1 

2(Az) 



JV-1 
J=2 



N—l 



Ho 



N 2 



(Ax) 2 



-2 1 
1 -2 



1 -2 1 
1 -2 



(100) 



(101) 



(102) 
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Figure 14: Heat equation: plots of Lyapunov functions 'HiAx (left) and "H.2^x (right) vs 
time, AVF integrator. 



and the obvious discretization of N\- With these choices, both discretizations 
yield identical semi-discrete equations of motion and therefore 'Hi and H2 are 
simultaneously Lyapunov functions of the semi-discrete system and therefore, 
the AVF integrator preserves both Lyapunov functions. 
Initial conditions and numerical data: 



x e [0, 1] 



N = 50, 



At = 0.0025. 



(103) 



Initial condition: u(x, 0) = x(l — x). 

This system is numerically illustrated in figure 14 where the monotonic 



decrease of the Lyapunov functions for the heat equation in ( 100 1 and ( 101 ) is 
shown. 



4. Concluding Remarks 

The concept of energy, i.e., its preservation or dissipation, has far reaching 
consequences in the physical sciences. Therefore many methods to preserve en- 
ergy, and several to preserve the correct dissipation of energy (e.g. |13U23j ). have 
been proposed for ordinary differential equations. Surprisingly, when partial dif- 
ferential equations are considered, a unified way to discuss the preservation or 
correct dissipation of energy is missing and similar ideas arc often developed 
from scratch (e.g. [TTJ [12]). In this paper, we have presented a systematic 
and unified way to discretize partial differential equations and to preserve their 
correct energy preservation, or dissipation, by the average vector-field method. 

For the equations treated in this paper, one can replace the average vector- 
field method by any energy-preserving B-series method, while retaining the ad- 
vantageous properties of energy preservation or dissipation. More generally, geo- 
metric integrators for Hamiltonian or non-Hamiltonian PDEs with non-constant 
matrix T> can be constructed using discrete gradient methods, cf. 23J. 
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